home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Power Programmierung
/
Power-Programmierung (Tewi)(1994).iso
/
magazine
/
progjour
/
1988
/
01
/
chilc
/
ssu.for
< prev
Wrap
Text File
|
1987-09-03
|
13KB
|
432 lines
C./ ADD NAME=SSU
C SUBROUTINE TO TEST ACCURACY OF ROOTS OF CUBIC EQUATION,
C
C SUBROUTINE ROOTST ( K, KAPPA, H )
C IMPLICIT REAL*8 (A-H,O-Z)
C REAL*8 K, K2, K3, KAPPA
C
C K2 = K*K
C K3 = K2*K
C T1 = DABS(2.*K3)
C T2 = H*K2
C DT = DABS(T1-DABS(KAPPA-T2))
C T = T1
C IF (DABS(KAPPA).GT.T) T = DABS(KAPPA)
C IF (DABS(T2).GT.T) T = DABS(T2)
C DT = DT/T
C IF (DT.LT..0001) RETURN
C PRINT 1000, KAPPA, H, K, DT
C1000 FORMAT (' ERROR IN ROOT, KAPPA = ',E12.4,' H = ',E12.4,'
C * K = ',E12.4,' DT = ',E12.4 )
C END
C
C SUBROUTINE TO CALCULATE THE COMPLETE ELLIPTIC INTEGRAL
C OF THE THIRD KIND,
C
SUBROUTINE CEL3( CE, KP, M, FP )
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 KC, M, M0, KP
KC = KP
P = FP
IF (KC*P.EQ.0.D0) GO TO 10
KC = DABS(KC)
E = KC
M0 = 1.D0
IF (P.LE.0.D0) GO TO 2
C = 1.D0
P = DSQRT(P)
D = 1.D0/P
GO TO 3
2 G = 1.D0 -P
F = KC*KC - P
P = DSQRT(F/G)
D = -M/(G*P)
C = 0.D0
3 F = C
C = C + D/P
G = E/P
D = 2.D0 * ( F*G + D )
P = G + P
G = M0
M0 = KC + M0
IF ( DABS(G-KC).LE.1.D-4*G ) GO TO 4
KC = 2.D0 * DSQRT(E)
E = KC * M0
GO TO 3
4 CE = 1.5707963D0 * (C*M0+D)/(M0*(M0+P))
RETURN
10 CE = (1.D18)**2
RETURN
END
C
SUBROUTINE SS( E, DE, S, DELTA, EC, BET )
C SIGNALS
C DLE = .0 AND DE NOT= 0, NORMAL,
C DLE > .0, COMPUTES COEFFICIENT OF LN(E-EL) NEAR E = EL,
C DLE < .0, COMPUTES COEFFICIENT OF (E-EL)**1/3 FOR E NEAR
C EL AND EL = EMID,
C DE = .0, GIVES END CORRECTIONS, S FOR E < OR E > EMAX,
C
IMPLICIT REAL*8 (A-H,O-Z)
real*8 darcos
DIMENSION S(3), X(3)
REAL*8 KAPPA, KAPRM, K(3), KL, NUM(3), K2, KC, KK, KC2
COMMON KAPPA, KAPRM, H, Q, RQ, RP, H2, R1
COMMON/A/NA, NS, KL, ISIG
COMMON/C/EPS, DLE
COMMON/D/S12, S11 /PI/PI
PI = DATAN(1.D0)*4.
DO 1 N = 1, 3
1 S(N) = 0.D0
DDE = .1*DE
E2 = EC - DELTA/3.D0
E3 = EC + DELTA/3.D0
EMIN = .0
IF (E2.LT.0.D0) EMIN = E2
EMAX = .0
IF (E3.GT.0.D0) EMAX = E3
EMID = E2
IF (EMID.EQ.EMIN) EMID = .0
IF (EMID.EQ.EMAX) EMID = E3
IF (DABS(E).LT.DDE) E = .0
IF (DABS(E-E2).LT.DDE) E = E2
IF (DABS(E-E3).LT.DDE) E = E3
B2 = BET**.66666666D0
B4 = B2**2
RP = B4 + 2.D0/B2
RQ = 18.D0/RP
IF (DLE.GT.0.) E = E-EPS
H = RP*E-2.D0*EC/B2
H2 = 6.D0*E - 4.D0*EC
Q = B4*(2.D0*H**2/3.D0)-BET**2*E**2-2.*((E-EC)**2+DELTA**2/9.D0)
Q = 3.D0*Q/(RP*B4)
R1 = 2.D0*H*(2.*BET**2+1.D0)/(RP*B2)
KAPPA = E*(E-E2)*(E-E3)
KAPRM = E*(E-E2) +E*(E-E3) +(E-E2)*(E-E3)
XP = KAPRM/RP
Y = ( Q + RQ*XP )
YP = ( R1 + 18.D0*H2/RP**2)
P23 = 2.D0*3.1415927D0/3.D0
IF (DLE.GT.0.D0) GO TO 50
C IF (DLE.EQ.0.D0) GO TO 7
C EL = EMID
C E = EL
C XA = ((EMID-EMIN)*(EMAX-EMID)*DE/2.D0)**.66666666D0
C A = DSQRT(3.D0)*XA
C RA = DSQRT(A)
C AK = DSQRT((1.D0-DSQRT(3.D0)/2.)/2.)
C AP = A+XA-XP
C CALL CEL1(KK, AK, IER)
C FJ1 = -KK/RA
C FJ2 = FJ1/(XP*RP)
C FJ3 = FJ2/(XP*RP)
C EK = .0
C GO TO 30
C
C NORMAL CASE
7 G = H/3.D0
IF (DE.EQ.0.) GO TO 40
IF (DABS(H**3).LT.DABS(KAPPA*1.E-10)) GO TO 4
C THREE ROOT REGION
C = 2.D0* KAPPA*(3.D0/H)**3 -1.D0
IF (C*C.GE.1.) GO TO 10
PHI = DARCOS(C)/3.D0
K(1) = DABS( G*(DCOS(PHI) -.5D0 ))
C
K(2) = DABS( G*(DCOS(PHI-P23) -.5D0 ))
C
K(3) = DABS( G*(DCOS(PHI+P23) -.5D0 ))
C
IF (ISIG.GT.1) PRINT 1000, (K(L), L = 1, 3 )
1000 FORMAT (' K S = ', 3E12.4 )
DO 3 L = 1,3
3 X(L) = K(L)**2
IF(DABS(E).LT.DDE.OR.DABS(E-E2).LT.DDE.OR.DABS(E-E3).LT.DDE)
* GO TO 18
K2 = (X(2)-X(1))/(X(3)-X(1))
AK = DSQRT(K2)
ALPHA2 = (X(2)-X(1))/(XP-X(1))
KC2 = (X(3)-X(2))/(X(3)-X(1))
KC = DSQRT(KC2)
FP = (XP-X(2))/(XP-X(1))
CALL CEL1( KK, AK, IER )
CALL CEL2( EK, AK, 1.D0, KC2, IER )
CALL CEL3( PIK, KC, K2, FP )
W = ( KAPPA - H*XP )**2 - 4.*XP**3
RTW = DSQRT(W)
WP = -2.*H2*(6.D0*XP**2-XP*H**2+H*KAPPA)/RP
RAC = DSQRT(X(3)-X(1))
X1P = X(1) -XP
X2P = X(2) -XP
X3P = X(3) -XP
FJ1 = -2.*KK/RAC
FJ2 = 2.*PIK/(RP*RAC*X1P) - PI/(RP*RTW)
FJ3 = -2.*RAC*EK/(W*RP**2)+.5D0*KK/(RAC*X1P*X2P*RP**2)
* -(6.D0*XP**2-XP*H**2+H*KAPPA)*PIK/(RAC*W*X1P*RP**2)
FJ3 = 2.*FJ3+PI*(6.D0*XP**2-XP*H**2+H*KAPPA)/(W*RTW*RP**2)
S(1) = YP*FJ2-H2*Y*FJ3
S(2) = 18.D0*FJ1/RP**2-Y*FJ2
S(3) = 2.*(X3P*KK/RAC-RAC*EK)*RP**2/3.D0
GO TO 20
C
C ONE ROOT REGION
4 K(3) = (DABS(KAPPA/2.))**.33333333D0
C
GO TO 11
10 TEMP = ( DABS(C) + DSQRT(C*C-1.D0))**.33333333D0
TEMP = DSIGN(TEMP,C)
K(3) = DABS( .5D0*( TEMP -1.D0 +1.D0/TEMP )*G)
C
11 K(1) = -1.
IF (ISIG.GT.1) PRINT 2000, K(3)
2000 FORMAT (' K(3) = ',E12.4 )
18 XA = K(3)**2
A2 = (6.D0*XA**2-XA*H**2+H*KAPPA)/2.D0
A = DSQRT(A2)
B1 = H**2/8. -XA/2.
K2 = ( A +B1 -XA )/(2.*A)
KC2 = (A -B1 +XA )/(2.*A)
KC = DSQRT(KC2)
AP = A+XA-XP
AM = A-XA+XP
C IF (DABS(AM/A).LT..001) AM = XP+XA*(DSQRT(3.D0+H/K(3))-1.)
XAP = XA-XP
ALPHA = AM/AP
FP = AP**2/(4.*A*XAP)
RT = DSQRT(K2+AM**2/(4.*A*XAP))
RA = DSQRT(A)
IF (K2.GT.0) GO TO 19
AK = .0D0
KK = PI/2.D0
EK = KK
PIK = KK/DSQRT(FP)
GO TO 17
19 AK = DSQRT(K2)
CALL CEL1 ( KK, AK, IER )
CALL CEL2 ( EK, AK, 1.D0, KC2, IER )
CALL CEL3 ( PIK, KC, K2, FP )
17 FJ1 = -KK/RA
AL2M = -(4.*A*XAP)/AP**2
A2OM = -AM**2/(4.*A*XAP)
AKK = KC2*ALPHA**2 +K2
IF (DABS(ALPHA).LT..01) GO TO 25
FJ2 = (AP*PIK/(2.*XAP) -KK)/(RA*AM*RP)
FJ3 = FJ2/(AM*RP) -(A*PIK/(2.*XAP**2) -4.*A**2*(A2OM*EK
* -AKK*KK/AL2M +(A2OM**2-K2)*PIK)/(ALPHA*AP**3*AKK))
* /(RA*AM*RP**2)
GO TO 30
25 FJ2 = KK/(2.*RA*RP*XAP)
IF (K2.GE..0001) TEMP = (KK-EK)/K2
IF (K2.LT..0001) TEMP = PI/4.
FJ3 = (.5D0*TEMP-KK)/(2.*RA*(RP*XAP)**2)
30 S(1) = YP*FJ2 -H2*Y*FJ3
S(2) = 18.D0*FJ1/RP**2 -Y*FJ2
S(3) = (AP*KK/RA -2.*RA*EK)*RP**2/3.D0
5 IF (DABS(E-EMID).GT.DDE.OR.DLE.LT.0.) GO TO 20
C
C CORRECT FOR DISCONTINUITY AT E = EMID
NUM(2) = -Q/KAPRM**2
NUM(1) = R1/KAPRM**2 +H2*NUM(2)/KAPRM
NUM(3) = RP/3.D0
G = 1.5707963D0*DABS(KAPRM)/(2.*K(3))
S11 = G*NUM(1)
S12 = G*NUM(2)
DO 8 L = 1, 3
8 S(L) = S(L) + G * NUM(L)
GO TO 20
C
C CALCULATE OUTSIDE S = S11, S12
40 RT = DSQRT((KAPPA-H*XP)**2-4.*XP**3)
S12 = -PI*Y/(RP*RT)
IF (ISIG.GT.3) PRINT 5000, E, S12
5000 FORMAT (' E = ',E12.4,' S12 = ',E13.6)
GO TO 20
C
C CALCULATE COEFFICIENT OF LN TERM,E NEAR EL
50 XC = (H/3.D0)**2
P = KAPRM - RP*XC
TEMP = -(DABS(P)/DSQRT(3.*XC))
C TEMP = TEMP*((2.+3.*EPS/DE)*ALOG(DABS(H*P*(DE+EPS)/(81.*XC**2)))
C * +(2.-3.*EPS/DE)*ALOG(DABS(H*F*(DE-EPS)/(81.*XC**2))) -6.)
P2 = P**2
S(3) = TEMP*RP/3.D0
S(2) = -TEMP*(Q+RQ*XC)/P2
S(1) = TEMP*R1/P2 +H2*S(2)/P
E = E + EPS
C
C COMMON TERMINATION
20 IF (ISIG.GT.3) PRINT 3000, E, ( S(L), L= 1,3 )
3000 FORMAT (' E = ',E12.4,' S S = ',3E14.8 )
C IF (IABS(ISIG).GT.1) PRINT 4000, FJ1, FJ2, FJ3
4000 FORMAT (' FJ S ',3E14.8 )
RETURN
END
C
REAL*8 FUNCTION DARCOS(X)
REAL*8 X,PI
COMMON/PI/PI
IF (X**2-.0001) 1, 1, 2
1 DARCOS = PI/2. - X - X**3/6.
RETURN
2 DARCOS = DATAN( DSQRT(1.D0-X**2)/X)
IF (X.LT.0.) DARCOS = DARCOS + PI
RETURN
END
C
C .............................................................................
C SUBROUTINE CEL1
C
C PURPOSE
C CALOULATE COMPLETE ELLIPTIC INTEGRAL OF FIRST KIND
C
C USAGE
C CALL CEL1(RES,AK,IER)
C
C DESCRIPTION OF PARAMETERS
C REC - RESULT VALUE
C AK - MODULUS (INPUT)
C IER - RESULTANT ERROR CODE WHERE
C IER=0 NO ERROR
C IER=1 AK NOT IN RENGE -1 TO +1
C
C REMARKS
C THE RESULT IS SET TO 1,E75 IF ABS(AK) GE 1
C THE RESULT IS SET TO 1,0E38 IF ABS(AK) GE 1
C FOR MODULUS AK AND COMFLEMENTARY MODULUS CK,
C EQUATION AK*AK*CK*CK=1,0 IS USED
C AK MUST BE IN THE RANGE -1 TO +1
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C NONE
C SPMAX
C
C METHOD
C DEFINITION
C CEL1(AK)=INTEGRAL(1/SQRT((1+T*T)*(1+(CKKT)**2)), SUMMED
C OVER T FROM 0 TO INFINITY),
C EQUIVALENT ARE THE DEFINITIONS
C CEL1(AK)=INTEGRAL(1/(COS(T)SQRT(1+(CK*TAN(T))**2)),SUMMED
C OVER T FROM 0 TO PI/2),
C CEL1(AK)=INTEGRAL(1/SQRT(1-(AK*SIN(T))**2),SUMMED OVER I
C FROM 0 TO PI/2), WHERE K=SQRT(1,-CK*CK),
C EVALUATION
C LANDENS TRANSFORMATION IS USED FOR CALCULATION,
C REFERENCE
C R,BULIRSCH, 'NUMERICAL CALOULATION OF ELLIPTIC INTEGRALS
C AND ELLIPTIC FUNCTIONS',HANDBOOK SERIES SPECIAL FUNCTIONS,
C NUMERISCHE MATHEMATIK VOL, 7, 1965, PP, 78-90,
C
C ............................................................................
C
SUBROUTINE CEL1(RES,AK,IER)
IMPLICIT REAL*8 (A-H,O-Z)
IER=0
ARI=2.
GEO=(0.5D0-AK)+0.5D0
GEO=GEO+GEO*AK
RES=0.5D0
IF(GEO)1,2,4
1 IER=1
C 2 RES=1.E75
2 RES = SPMAX(1)
RETURN
3 GEO=GEO*AARI
4 GEO=DSQRT(GEO)
GEO=GEO+GEO
AARI=ARI
ARI=ARI+GEO
RES=RES+RES
IF(GEO/AARI-0.9999)3,5,5
5 RES=RES/ARI*6.283185D0
RETURN
END
C
C ..........................................................................
C
C SUBROUTINE CEL2
C
C PURPOSE
C COMPUTE THE GENERALIZED COMPLETE ELLIPTIC INTEGRAL OF
C SECOND KIND,
C
C USAGE
C CALL CEL2(RES,AK,A,B,IER)
C DESCRIPTION OF PARAMETERS
C REC - RESULT VALUE
C AK - MODULUS (INRUT)
C A - SONSTANT TERM IN NUMERATOR
C B - FACTOR OF QUADRATIC TERM IN NUMERATOR
C IER - RESULTANT ERROR CODE WHERE
C IER=0 NO ERROR
C IER=1 AK NOT IN RENGE -1 TO +1
C
C REMARKS
C FOR ABS(AK) GE 1 THE RESULT IS SET TO 1,E75 IF B IS
C POSITIVE, TO -1,E75 IF B IS NEGATIVE,
C FOR ABS(AK) GE 1 THE RESULT IS SET TO 1,0E38 IF B IS
C POSITIVE, TO -1,0E38 IF B IS NEGATIVE,
C SPECIAL CASES ARE
C K(K) OBTAINED WITH A = 1, B = 1
C E(K) OBTAINED WITH A = 1, B = CK*CK WHERE CK IS
C COMPLEMENTARY MODULUS,
C B(K) OBTAINED WITH A = 1, B = 0
C D(K) OBTAINED WITH A = 0, B = 1
C WHERE K, E, B, D DEFINE SPECIAL CASES OF THE GENERALIZED
C COMPLETE ELLIPTIC INTEGRAL OF SEOOND KIND IN THE USUAL
C NOTATION, AND THE ARGUMENT K OF THESE FUNCTIONS MEANS
C THE MODULS,
C
C SUBROUTINE AND FUNCTION SUBPROGRAMS REQUIRED
C NONE
C SPMAX
C
C METHOD
C DEFINITION
C RES=INTEGRAL((A+B*T*T)/(SQRT((1+T*T)*(1+(CK*T)**2))*(1+T*T))
C SUMMED OVER T FROM 0 TO INFINITY),
C EVALUATION
C LANDENS TRANSFORMATION IS USED FOR CALCULATION,
C REFERENCE
C R,BULIRCH, 'NUMERICAL CALCULATION OF ELLIPTIC INTEGRALS
C AND ELLIPTIC FUNCTIONS' , HANDBOOK SERIES SPECIAL FUNCTIONS,
C NUMERISOHE MATHEMATIK VOL, 7, 1965, PP, 78-90,
C
C ..........................................................................
C
SUBROUTINE CEL2(RES,AK,A,B,IER)
IMPLICIT REAL*8 (A-H,O-Z)
IER=0
ARI=2.
GEO=(0.5D0-AK)+0.5D0
GEO=GEO*AK + GEO
RES=A
A1=A+B
B0=B+B
IF(GEO)1,2,6
1 IER=1
2 IF(B)3,8,4
C 3 RES=-1.E75
3 RES=-SPMAX(1)
RETURN
C 4 RES=1.E75
4 RES=SPMAX(1)
RETURN
5 GEO=GEO*AARI
6 GEO=DSQRT(GEO)
GEO=GEO+GEO
AARI=ARI
ARI=ARI+GEO
B0=B0+RES*GEO
RES=A1
B0=B0+B0
A1=B0/ARI+A1
IF(GEO/AARI-0.9999)5,7,7
7 RES=A1/ARI
RES=RES+0.5707963D0*RES
8 RETURN
END
REAL*8 FUNCTION SPMAX(K)
SPMAX=1.0D38
RETURN
END